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' Abstract: We present a simple method to automatically evaluate arbitrary tree-level am- 

■ plitudes involving the production or decay of a heavy quark pair QQ in a generic 2S+i^[^>8] 

. state, i.e., the short distance coefficients appearing in the NRQCD factorization formal- 

ism. Our approach is based on extracting the relevant contributions from the open heavy 
^ , quark-antiquark amplitudes through an expansion with respect to the quark-antiquark rel- 

H ' ative momentum and the application of suitable color and spin projectors. To illustrate 

the capabilities of the method and its implementation in MadGraph a few applications to 
quarkonium collider phenomenology are presented. 
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1. Introduction 

Bound states of heavy quarks, such as J/i^, T or are unique systems where our under- 
standing of QCD both at the perturbative and non-perturbative level can be tested. Since 
the mid seventies, when J/ip was discovered and asymptotic freedom invoked to explain its 
narrow width, an impressive amount of experimental data have provided continuous chal- 
lenges to the theoretical community. To date several open questions remain to be answered 
that concern both production and decay mechanisms Q. 

One of the outstanding issues is to what extent a non-relativistic treatment of a quarko- 
nium system is a useful starting point to describe the production and decay of bound states 
of heavy quarks. This is best tested in the case of charmonium, where the mass is not too 
far from Aqcd- The appeal of such an approach, is the possibility of embedding it in a 
rigorous framework, non-relativistic QCD (NRQCD) |Q, that allows one to make consistent 
theoretical predictions which can systematically be improved. NRQCD is an effective field 
theory that can be matched to full QCD. Effects of order > rriQ are encoded in short 
distance coefficients that can be estimated by using perturbation theory. Low < rriQ 
effects, typically those involving hadronization, are factorized into long-distance matrix 
elements, which are universal, can be organized hierarchically in powers of v and can be 
either measured on lattice simulations or extracted from experimental data. 
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Despite its theoretical appeal and undeniable successes, not all of the predictions of 
the NRQCD factorization approach have been firmly established. Consider for example 
the color-octet contributions in the J/ip production: they seem to play a dominant role 
at the Tevatron and in 77 collisions, but they look marginal at e^e~ at low energy, in 
photoproduction at HERA Q and in fixed-target experiments [Q]. Another example is the 
of J/^jj 's at high-pr at the Tevatron. NRQCD predicts a sizeable transverse polarization 
in contrast to the latest data that now clearly indicate that J/tp^s are not transversely 
polarized Q. 

Given such a puzzling scenario, it is mandatory to re-examine in details the key ob- 
servables and the corresponding available theoretical predictions and try to systematically 
improve on them. Two main directions can be followed. 

The first is to calculate predictions for cross sections and relevant distributions at 
the next-to-leading order (NLO) accuracy in the strong coupling expansion. Even for a 
small number of final state partons these calculations are technically very challenging and 
only a few are currently available. For example, there are only three calculations relevant 
for describing the pT of the quarkonium state at NLO at colliders, i.e., those whose born 
process is a 2 ^ 2: the 77 + 7 + X Q, the pioneering — > ^5^^^ + X and 

the recently evaluated pp "^S^^ + X [^. In e^e~ collisions the cross section to produce 

cc is the only inclusive charmonium cross section known at NLO 

The second way is to extend the set of the typical observables analyzed by studying 
new processes and also predicting more exclusive observables. To this aim, we propose 
to use a multi-purpose matrix element based generator. This kind of codes, which have 
become widely available recently, have dramatically reduced the burden of performing 
phenomonological studies at colliders. They automatically create the matrix elements cor- 
responding to a given process and then generate unweighted parton-level events that can 
be passed to standard Monte Carlo programs such as Pythia or Herwig for show- 
ering and hadronization, and eventually to detector simulation. The advantages of such 
an approach are numerous. The most important one is that any theoretical description 
(model) can be directly used in the experimental analyses and any kind of realistic observ- 
able can then be studied both at the theoretical and experimental (detector) level, within 
the same framework. Second, since the computation of the matrix element is automatic, 
virtually any process can be simulated at tree level accuracy, opening up the possibility of 
straightforward experimental analyses of potentially interesting ideas. 

The purpose of this work is to present the first step in this direction, i.e., an algorithm 
for the automatic matrix element generation of amplitudes involving the production (or 
the decay) of a short-distance quarkonium state with arbitrary spin, angular momentum, 
and color quantum numbers. The algorithm has been implemented in the matrix element 
generator MadGraph |12|. Inclusion of the matrix elements into the automatic phase 
space integrator and event generator MadEvent |13, 14] is ongoing and will be presented 
in a second publication. The plan of the paper is as follows. In Section || we describe 
the approach followed to project open quark-antiquark amplitudes into those for a given 
^■^^^Lj'^' state and the validation performed. In Section ^ we present three applications 
to collider phenomenology. We draw our conclusions in the final section. 
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2. Quarkonium amplitudes in NRQCD 



In the framework of Non Relativistic QCD, the cross section associated to a heavy quarko- 
nium production process can be factorized as 

criij ^Q + X) = Y, Hij - QQin) + X)(OS(n)) , (2.1) 

n 

where 

• a{ij — > QQ{n) + X) is the short distance cross section, 

• {0^{n)) is the long distance matrix element. 

The short distance cross section a{ij QQ{n) + X) is related to the creation of a heavy 
quark pair in a quantum state n and it can be computed within perturbation theory, 
by expanding the amplitude in powers of as and v, i.e., the strong coupling constant 
and the relative velocity v. The sum over n in Eq. (|2.1| ) a priori runs over all possible 
non-pertubative transitions QQ{n) Q and goes therefore beyond the color-singlet ap- 
proximation, where only the perturbative state with the same quantum numbers as the 
physical state is included. In particular, the intermediate QQ state can be in a color-octet, 
the color being neutralized at later times. 

In practice, any quarkonium squared amplitude is obtained by means of an expansion 
in powers of the relative velocity v between the heavy quarks. At leading order in v, the 
intermediate QQ state can be specified by the spectroscopic notation 

n= 2^^+14^], (2.2) 

where S identifies the spin state of the heavy quark pair, L the orbital momentum state, 
J the total angular momentum state, and c = 1, 8 the color state. 

To evaluate the short distance coefficients we use the general algorithms for Feynman 



diagrams based computations available in MadGraph |12]. MadGraph can provide partonic 
helicity amplitudes for any process in the Standard Model and beyond at the tree level. 
In the case of quarkonium, we start from the helicity amplitude information for the open 
quark-antiquark production and then apply several projection subroutines that select the 
specific quantum state numbers, ^^~^^L^j\ The projectors are universal and depend only 
on the quantum numbers to be selected and not on the specific process requested. In fact, 
only the color projection has to be performed within MadGraph itself since it changes the 
way the amplitudes are calculated (and in particular selects which diagrams contribute) 
while spin projections and velocity expansions are performed at runtime. Our method 
entails that any squared matrix amplitude — > QQ(^^~^^ L^j^) + X)\'^ can be generated 
in any of the models available in MadGraph. Limitations are typical to Feynman-diagram 
based matrix element generators and are connected to the factorial growth in the number 
of diagrams. Let us stress a few important points. First, our implementation also allows 
the calculation of the relativistic corrections for S-wave state production. Second, the 
employed projectors are general enough to handle the case where the quark and the anti- 
quark are of different fiavor, and therefore have a different mass. In our approach spin 
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correlation between the quarkonium angular momentum and its decay products are fully 
included and therefore the information on the quarkonium polarization is kept. In Section 



2^ , we discuss two methods to disantangle the longitudinal and transverse parts of the 
cross section for the transition QQ{^S^^^) Q{J^^ = 1""). 

2.1 Projection method and its implementation 

In this section, we discuss how to combine the open-quark helicity amplitudes generated by 
MadGraph in order to select a specific configuration ^s+i^W ^-^^ heavy-quark pair. Our 



method can be seen as a generalization and automatization of that proposed in Refs. |15 



First we address the computation of the color factors. In MadGraph, the color struc- 
ture of the parton-level amplitude is obtained by organizing the amplitude into gauge 
invariant subsets, corresponding to different color flows in the large- A'^c limit |17]. The 
total amplitude is then expressed as a sum over the color flows 

M = Y,fm{aA---)Am, (2.3) 

m 

where the /m's are color factors (orthogonal in the large Nc limit) that depend on color 
indices a,b, . . . carried by gluons and quarks, and the A^s are gauge invariant quantities, 
the so-called partial amplitudes. These amplitudes are computed in MadGraph with the 
help of HEL AS subroutines ||l8| . The color factors are evaluated at the level of the squared 
amplitude, each color index being contracted: 



\M\^ = {Al,A*,,...) 



/l/i fU2 ... I A2 \ . (2.4) 




In MadGraph, the decomposition of the partonic amplitude into partial gauge-invariant 
amplitudes, as well as the evaluation of the matrix of color factors, are computed automat- 
ically for any process at tree-level. 

In order to compute the color structure related to the process ij —>■ QQ(^^~^^L^j) + X, 
the quark and anti-quark are required to be in a given color state [c] . Given that 3 3 = 
1 © 8 the pair can be either in a singlet or octet state that can be obtained by using the 
projectors [^] 

where i and j are the color indices of the heavy quarks. These projectors have been imple- 
mented into MadGraph. The color structure for the quarkonium production amplitude is 
then computed using the same technique as for the open quark production. 

Next we consider the projection of the amplitude onto a definite total spin state of 



Ci = ^, Cs = V2iT%„ (2.5) 



the heavy quark pair. The spin projectors were first derived in Ref. |20|, and we employ 



the normalization of Ref. [O] . Since the code is also designed for the production of heavy 
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quarkonia of mixed flavours, we distinguish the mass of the quark (rn-i) from that of the 
anti-quark {1712)- pi {p2)i the momentum of the heavy quark (anti-quark) can be expressed 
in terms of the total momentum P and the relative momentum p: 

mi 7712 

Pi = ; P + P, P2 = \ P-p. (2.6) 

777-1 + "1-2 "T-l + "^2 

The partial helicity amplitudes can be combined to yield a specific spin state for the heavy 
quark pair: 

M{ij ^QQ{S,\)+X)= Y,N{X\Xi,\2)M{ij^Q{Xi)Q{\2)+X), (2.7) 

where the Clebsch-Gordant coefficients A^(A|Ai,A2) can be written in terms of the heavy 
quark spinors ^, 

iV(A|Ai, A2) = Ji v{p2, \2)Tsu{pi, Ai) , (2.8) 

V8mim2 

with Ts=i = 6^7^^ for a spin-one QQ state (e^ being the polarization vector), and Ts=o = 75 
for a spin-zero QQ state. Note that this expression is accurate only at leading order in the 
relative momentum p. The helicity amplitudes M.{ij — > Q{Xi)Q{X2) + X) are computed 
by MadGraph in a specific Dirac representation employed in the HELAS subroutines |l8| . 
These subroutines can also be invoked for the evaluation of the current in the rhs of 



Eq. (2.8). The spin projection formula in Eq. ( p.TD can then be easily incorporated in a 



generic numerical algorithm. 

We then consider the projection onto a specific orbital state. For 5-wave state produc- 
tion at leading order in one can simply set p = in the short distance amplitude. For 
P-wave state, the leading order contribution in v is given by the derivative of the amplitude 
with respect to the relative momentum. We approach this derivative numerically by the 
quotient 

M{/\p')-M{Q) 

in the quarkonium rest frame^. 

Eventually, for spin-one P-wave states, the spin index is combined with the orbital 
index in order to select a given angular momentum state, J = 0, 1 or 2. 

2.2 Relativistic correction to S'-wave state production 

To evaluate the relativistic corrections to S'-wave state production, the relative momentum 
p must be kept different from zero in the amplitude. In this case, we need an expression 



of the spin projectors accurate to all orders in p. According to the result of Ref. |21|, the 
Clebsch-Gordan coefficients can be written as 

7V(A|Ai,A2) = -i;(p2,A2)^4^r5n(pi,Ai) , (2.10) 

\/2(P + m) 4±!/ 



^We use the normalization u{pi, \i)u{pi, Ai) = 2mi = — u(pi, \i)v{pi, A;). 

^Wherever we consider a non-zero relative momentum, we replace by Et (the energy of the quark i) 



in Eq. (2.6) in order to enforce the mass-shell condition for the heavy quarks. 
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where E = In this case we assume that the quark and the anti-quark are of the 

same flavour (mi = m2 = m). The contribution to an 5"- wave conflguration is selected by 
projecting the amphtude onto the spherical harmonic Yj^"^^ . This amounts to averaging 
the amplitude over the direction of the relative momentum. The short distance amplitude 

2 

is then expanded in powers of each term in the expansion being associated with an 
independent non-perturbative long distance matrix element. 

It has been shown recently that the relativistic corrections arising from the wave func- 



tion for S'-wave state can be resummed to all orders inv |22, Thanks to the following 
relation among the long distance matrix elements 



(0|xtK(-v2)"v|Q) _ /(0|x^k(-v2)V'|q) 



(oixt^^VIQ) V (oix^AcVIQ) 



(2.11) 



where k is either 1 for spin-zero QQ states or the Pauli matrix a for spin-one QQ states, 
the cross section can be expressed in term of one long distance matrix element: 



a 



>mx^^^m){Q\^^^^x^ , (2.12) 



where the short distance coefficient has been evaluated at the relative momentum 

The inclusion of the relativistic corrections arising from the wave function has been imple- 



mented in the code, using Eq. ([2.121 ). 



2.3 Polarization of vector-like quarkonium states 

The polarization of vector-like quarkonium states is an interesting exclusive quantity that 
can provide further information on the production mechanism since it can be measured 
accurately at colliders. On the theoretical side, the transverse and longitudinal parts of 
the cross section can be disentangled by using the explicit basis of polarization vectors 
(e^, en, £7^2)5 defined by the equations 

Ei.Sj = -Sij, Si.P = 0, elP = 0. (2.14) 

This representation is obviously Lorentz-frame dependent. Usually, the three- vector P is 



evaluated in the lab frame in Eq. (2.14) but other useful frames exist |24|. 

At hadron colliders a ^Si quarkonium state is observed through its decay into leptons. 
The polarization of the quarkonium can indeed be determined by analyzing the angular 
distribution of the leptons. Defining as the angle between the direction in the quarko- 
nium rest frame and the quarkonium direction in the laboratory frame, the normalized 
angular distribution /(cos 6) is 

I{cos9)= ^ (1 + acos^g), (2.15) 
2[a + 6) 

where the relation between a and the polarization state of the quarkonium is 

a = ^^^^ . (2.16) 
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In fact, a more general way to correctly account for the polarization information is to 
decay the vector-like quarkonium state into leptons at the matrix element level. This can 
be achieved by replacing the polarization vector of the quarkonium by the leptonic current 

Ui-{ki,Xi)^^Vi+{k2,X2) ■ (2.17) 



8■m^/^T 

and has been set up as an option in our code. Once the vector-like quarkonium decay 
is required by the user, events are generated that include the information on the lepton 
momenta too. This method has two main advantages. First it is directly connected to 
the experimental analysis. The acceptance of the detector may affect the relation between 
the measured value of the parameter a and the perturbative polarized cross sections ajn, 
gl. However, even when cuts are applied, the angular distribution of the leptons is the 
physical observable that can be used to compare theory and measurements. Second, it is 
Lorentz invariant and therefore does not impose an a priori choice for the frame where the 
polarization information is extracted. 

2.4 Validation and examples 
2.4.1 Preliminary checks 

In order to gain confidence in our implementation, we have performed several checks listed 
below. 

• Gauge invariance has been systematically verified for all the processes considered by 
using longitudinal polarization for the external photons or gluons. 

• Numerical cancellation among diagrams has been checked for amplitudes that vanish 
by symmetry considerations: 

A{^S^Q^ + {2k + 1) photons) = 0, 
^(^5^ + 2A; photons) = 0, 
^(ipji] + 2/c photons) = 0, 
^(3p|i] _^ 2k photons) = 0, 
A{^P^^l + {2k + 1) photons) = 0, 

with k = 1,2,3. 

• We have compared our numerical amplitudes against the analytic results point-by- 
point in phase space for {i,j, k = quarks or gluons): 

1. ij Qk for all S- and P-wave states, both color-singlet and color-octet ampli- 
tudes, i20lll2ll2f 



2. ij — > QV, with V = Z,W for the relevant S- and P-wave states, both color- 
singlet and color-octet amplitudes p7| ]; 

3. ij — > Q4>, 4> being a scalar or pseudo-scalar for the relevant S- and P-wave 
states, both color-singlet and color-octet amplitudes |28|; 
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ifb) 






ip[i] 


3p[l] 


3p[l] 


3p[l] 




1.69 • 10^ 


2.37 • 10^ 


1.78 • 102 


1.08 • 102 


2.23 • 102 


2.40 • 102 


77@LEP II 


0.519 


5.23 


0.162 


2.69 • 10-2 


5.80 • 10-2 


0.266 


7P@HERA 


3.66 • 10^ 


1.76 ■ 10^ 


85.4 


21.7 


51.7 


2.02 • 102 


pp@LHC 

qq 


3.94 • 10"^ 


9.83 • lO'^ 


5.20 • 10^ 


1.82 • 10^ 


4.40 • 10^ 


1.05 • lO'^ 


1.37- 10^ 


8.34 • 10^ 


2.95 • lO'' 


1.10 • 10^ 


2.69 • 10^ 


7.26 • 10^ 


pp@Tev II 


2.56 ■ 10^ 
2.64 • 10^ 


6.30 • 10^ 
1.62 • 10^ 


3.28 ■ 10^ 
5.70 • 10^ 


1.24 ■ 10^ 
2.12 • 10^ 


2.83 • 10^ 
5.21 • 10^ 


6.62 ■ 10^ 
1.41 • 10^ 



Table 1: Cross sections (femtobarn) for the color-singlet contributions to the inclusive production 
of Be- 



ifb) 


1g[8] 

Dq 


35.J8] 


lp[8] 


3p[8] 


3p[8] 


3p[8] 


e^e~@mz 


1.59 


2.22 


0.167 


0.102 


0.209 


0.225 


77@LEP II 


4.87- 10"^ 


4.90 • 10-3 


1.52 • 10-"^ 


2.52 • 10^^ 


5.44 • 10~^ 


2.49 • 10-^ 


759@HERA 


1.18 


8.46 


0.506 


7.64 • 10-2 


0.244 


1.61 


pp@LRC 

qq 


4.11 • 10^ 


1.79 • 10^ 


1.17 • 10^ 


1.38 • 10^ 


6.23 • 10^ 


2.29 • 10^ 


1.03 • 10^ 


7.03 • 10^ 


251 


30.7 


174 


616 


pp@TeY II 

qq 


2.85 • 10^ 
199 


1.26 • 10^ 
1.37 • 10^ 


8.13 • 10^ 
48.6 


9.80 • 102 
5.95 


4.25 • 10^ 
33.7 


1.57- 10^ 
120 



Table 2: Cross sections (femtobarn) for the color-octet contributions to the inclusive production 
of Be- 

The agreement for S wave amplitudes is at the machine precision, while for P-waves, 
which are obtained through a numerical derivative, it is typically at the 10-^ level. 

2.4.2 Example 1: Be production at various colliders 

Be production has been extensively studied at various colliders and the many available 
results provide us with a further testing ground for our code. 

We present the results for inclusive cross sections at several colliders for all the S and 
P, both color-singlet and color-octet, states. Since our purpose is to provide reference 
numbers, we use a very simple set of input parameters common to all processes: 

• mb = 4.9 GeV, rric = 1.5 GeV, ruB, = 6.4 GeV 

• (O (^^+^4^')) = (2J + 1) 0.736 GeV^ 

• {O = (2J + 1) 0.287 GeV^ 

• (0[81) = O.Ol(OW) 

• Mepa = IJ-F = IJ'R= 12.8 GeV 
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(a) 



(b) 



Figure 1: Higgs decays into T: in association with a bb pair (a) and with a photon (b). 



• asM = 0.189 



• pdf set: cteq611 

• aEM = 1^ 

Our results are summarized in Tables 1 and 2. For the processes gg Bcbc we 



However, some of the results in the table are new. 
2.4.3 Example 2: Higgs decay into T + X. 

In this section we illustrate through a very simple example how a multi-purpose matrix 
element code allows to quickly assess the phenomenological relevance of new ideas. 

A light Higgs of iriH ^ 140 GeV in the Standard Model as well as in extensions such as 
in SUSY or in generic 2HDM, decays predominantly into hb. Assuming a 2 — > 1 production 
mechanism at hadron colliders, such gg or bb fusion, a two 6-jet signature is obtained. Such 
events are extremely difficult to trigger on and to select from the enormous QCD two-jet 
background. In fact, in this mass range, the golden channel for discovery is H ^ which 
has a branching ratio of the order of only 10~^ in the SM but has a very clean signature: 
it can be triggered on and the invariant mass determination is very accurate, giving the 
possibility of reconstructing a small BW peak over a very large background. However, in 
some BSM scenarios, such as for example at large tan/? in SUSY or in a generic 2HDM, 
the branching ratio to two photons is highly suppressed and other decay modes need to be 
considered. 

The presence of a T in the final state could help both in the triggering (through the 
/i+yU^ pair) and the invariant mass reconstruction and give a viable discovery mode. Let 
us consider two cases. As the simplest decay mode, consider H — > Tbb, Fig |^ (left). This 
is analogous to Z — > Tbb, which can also be described by the fragmentation function 
approach. The exact matrix element calculation gives: 



compared with Refs. p9| , |30| , and found agreement for all intermediate states. For color- 
singlet 5- wave state production in jg interaction, our results agree with those of Ref. pT|| . 



r{H Tbb) 



= 1.56 • 10 



(2.18) 
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(a) 



(6) 



7 WVW 




wwv 




(^) 

Figure 2: J/i/j direct pliotoproduction. Production mechanism (a) proceeds via a ^sf^ state, while 
the photon associated (b), the cc associated and the multi-gluon (d) mechanisms can proceed via a 
color-singlet transition. 

where we used mn = 120 GeV, (Ot(^5'P)> = 9.28 GeV^ and as = 0.118. 

Another, potentially cleaner mode, which could allow for an excellent invariant mass 
resolution, is the decay H — > 7T, Fig || (right) Note that the crossed decay T — > jH has 
been considered before in the literature |33, 



search mode for a light scalar 



or pseudoscalar and it is known at the NLO accuracy |3g]. A straightforward calculation 



gives: 



T{H T7) 
T{H ^ 66) 



6.17-10" 



(2.19) 



where we used mn = 120 GeV, {Or{^S^^^)) = 9.28 GeV^ and oem = 1/137. 

The BR's for both decay modes in the SM (and also in SUSY, since possible large 
tan/3 enhacements cancel in the ratio of the widths), are found to be rather small, also 
considering that the branching ratio of the quarkonium state into leptons has yet to be 
included. 



3. Applications 

3.1 77 J/ij + X at LEP II 



J/ip production in 77 collisions at LEP II (\/5 = 196 GeV) has been studied in 
Ref In that paper, the direct (77), the single resolved (7(7) and the double resolved 
{gg) components have been taken into account. However, the only direct contribution that 
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Figure 3: Transverse momentum distributions for J/ip production in 77 collisions 



has been considered is the color-octet transition 

77-cc(35pl)5. (3.1) 

We now address the evaluation of the color-singlet J/ip production in (unresolved) 77 
interactions. Due to color conservation, the final states cc ^^S"!^^^ g is forbidden. But 
if we turn the gluon into a photon, the corresponding a|;j^,^-process is allowed. Charge 
conjugation conservation forbids the final state cc ^^'S'[^^^ gg- So the only color-singlet 
process at o'^m'^'s associated production. Finally, the channel J/ip+ light partons 

occurs at order Q;|;^a|. So the Born- level contributions for the color-singlet production 
are given by 

77 - cc 7, 77 ^ cc cc, 77 - cc ggg. (3.2) 

Note also that the final state cc ^^"S^^'^ gqq is also allowed, but is not finite: the fragmen- 
tation part 7 — > gg — > 5 is effectively included in the single-resolved contributions. 

Using the Weiszacker- Williams approximation at LEP2, as = 0.189, rrij^^ = 2mc = 3 
GeV, OEM = 1/137, (O (^^W^) = i.ig QeV^, {O (^sf )) = 1.06 • lO'^ GeV^, we obtain 
the Pt spectra displayed in Fig. ^ for the color-singlet and color-octet processes. We find 
that the associated cc production is by far the dominant one for Pt{J /^) > 5 GeV. 

3.2 e+e" -^ric + X at ^= 10.6 GeV 

The production of an r]c state from e^e~ collisions has been studied in the context of 
associated production [^]. Using the inputs as = 0.26, aEM = jgyj {O (^^^q^^^) = 0.387 
GeV^, we obtain 

(7(e+e- ^ 7/c + cc) = 58.7 fb, (3.3) 
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Figure 4: Representative diagrams for rjc electroproduction via color-singlet transition: at order 
Og in association with a cc pair (a) , and at order with three gluons (b) and with a light quark 
pair (c). 

in agreement with the result obtained in Ref. |]38|| . 

The contribution to rj^ production coming from the associated production with light 
partons has not been studied so far. The reason is that for color-singlet production, ne- 
glecting the exchange of an off-shell Z, there is no such processes at order (x^m'^'s^ '^^^ 
to color and charge conjugation conservation. However, at order a^j^ag the following 
processes can occur: 

e+e" Vcggg, e+e" rjcqqg , (3.4) 

which are finite and could give a non-negligible cross section due to the larger phase space 
available. This calculation is straightforward. The code generates the matrix elements for 
the subprocesses in ( |3.4D , whereas an analytical computation would have been cumbersome. 
Using a very simple phase-space generator, we checked that the channel r]c+ light partons 
constitutes approximately 10% of the inclusive cross section. The plot in Fig. ^ displays 

21 f I 

the differential cross section with respect to z = the fraction of momentum taken 

by the rye- The total cross sections are given by 

cr(e+e" r]cggg) = 3.72 fb, cr(e+e~ rjcqqg) = 1.63 fb (3.5) 

As a result, the associated production is dominating the r](. production 

a{e^e~ — > rjccc) 



91.6% (3.6) 



It is interesting to note that the corresponding predictions for J/ip production are 
quite different. A similar calculation gives: 

C7(e+e~ ^ J/V'cc) = 148 fb, a{e+e- ^ J /il^gg) = 266 fb , (3.7) 

and therefore, at leading order, 

a{e'^e~ — > J/ipcc) 



a{e+e- J/^X) 



35.7% . (3.8) 
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In contrast to the above LO estimate, data indicate that J/tj: is produced in association 
with charm most of the time 11391: 



4t EXP = 82 ±0.15 ±0.14%. 3.9 

It would be interesting to have the analogous measurement for the r/c too. 
3.3 ^ T ± X at the Tevatron 

As a last apphcation of our code, we consider the contributions coming from higher order 
processes to the inclusive production of a T(^S'[^^) at the Tevatron. In particular, we focus 
on the prediction for pT distribution. 

It is easy to verify that the partonic differential cross section at LO, Fig. ^(a), falls off 

very quickly with the pT, as l/pj', due to the presence of two very off-shell heavy quark 

„ [11 

propagators in all contributing diagrams. At NLO S\ production in association with 
gluons features diagrams with only one off-shell heavy quark whose contribution scale as 
1/pjn, Fig. P(b). At NNLO both fragmentation contributions. Fig. |^(c), and high-energy 
enhanced contributions, Fig. ^(d), appear. While a full NNLO calculation including vir- 
tual contributions is beyond our present technical capabilities, one can argue that the 
terms with a different scaling, such as l/pj^ and 1/p^ appear at the first time at tree-level 
at order Og and a| respectively, and therefore they are finite and are not affected by the 
(missing) virtual contributions taking place at the same order. ^ Even at tree-level, the 
calculation of multi-parton amplitudes is very challenging and an automatized computa- 
tional method is mandatory. For instance, the amplitude for gg — > '^5'[^' ± 2>g involves 
thousands of Feynman diagrams which reduce to several hundreds when the color and spin 



^This approach can be explicitly validated for the \/pt terms using the full NLO calculation [g]. A 
detailed phenomenological analysis is in progress [M. 
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Figure 6: Tree- level contributions to T hadroproduction via a color-singlet transition (a) at LO 
a|, (b) at NLO a|, (c,d) at NNLO a|. 

projections are applied. Another complication comes from the large number of different 
parton-parton subprocesses, whose bookkeeping becomes quite cumbersome for large mul- 
tiplicities. The above difficulties are easily solved by MadGraph. The identification of the 
relevant partonparton — > '^s'^'^ + Spartons subprocesses and the generation of the corre- 
sponding amplitudes is matter of a few seconds. In Fig. M we show the different scaling 

„ [11 

in pt for the various gluonic amplitudes gg — > S[ + ng, with n = 1,2,3. A minimal 
invariant mass for any gluon (initial or final state) pairs, sij > 4m^, is required to avoid 
the phase space regions where the matrix element is singular. The pT shapes, however, 
depends weakly on the cutoff. 

4. Conclusions and Outlook 

Quarkonium physics is a rich and active field of research, where our understanding of 
QCD can be challenged and tested. Collider data, in particular, suggests that our present 
description of the quarkonium production mechanism is not fully satisfactory yet. Given 
the importance that T and J/il^ production will have at the LHC not only as calibration 
tools, but also for bottom and even top physics, it is mandatory to improve the reliability 
and the flexibility of our predictions. In this work we have completed the first step towards 
the development of a multi-purpose Monte Carlo generator for quarkonium physics, i.e., a 
code for the automatic evaluation of any tree-level matrix element involving a heavy-quark 
pair in a definite spin and color state. To illustrate its possibilities we have presented 
several applications to collider physics phenomenology. 

The following step will be to interface the matrix element generator to MadEvent for 
the automatic phase space integration and event generation. With such a code, events at 
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Figure 7: Transverse momentum distribution for T production at the Tevatron, Run II. A cut on 
the minimal invariant mass of each pair of initial and final state gluons, Sij > 4m^, is applied to 
avoid soft and coUinear divergences. Only amplitudes involving gluons are included. Factorization 
and renormalization scales are fixed to fip = /i/? = 2™?,, with mi, = 4.75 GeV, {0{Y)) = 9.28 GeV'^. 
The PDF set is CTEQ6M and as(2mb) = 0.180. 



the parton-level for any process of interest can be generated and then passed through Pythia 
or Herwig for the parton shower and hadronization and finally to detector simulation. This 
will allow not only a much wider range of possible studies, such as for example pattern of 
extra radiation in quarkonium events, but also more direct flow of information from theory 
to experiment in testing new ideas or new approaches (such as different scaling rules for 
the various non-perturbative matrix elements). Work in this direction is in progress. 
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